Towards Better Integrators for Dissipative Particle Dynamics Simulations 



O 

o 
o 

(N 
-i— > 

o 
O 



Gerhard BesokpB, Ilpo Vattulainen 1 , Mikko Karttunen 2 , and James M. Poison 3 
1 Department of Chemistry, Technical University of Denmark, Building 207, DK-2800 Lyngby, Denmark 
2 Max Planck Institute for Polymer Research, Theory Group, PO Box 314-8, D-55021 Mainz, Germany 
3 Department of Physics, University of Prince Edward Island, 550 University Avenue, Charlottetown, PEI, Canada CIA £P3 
(Accepted for publication in PHYSICAL REVIEW E as Rapid Communication — tentative publication issue: 01 Dec 2000) 

Coarse-grained models that preserve hydrodynamics provide a natural approach to study collective 
properties of soft-matter systems. Here, we demonstrate that commonly used integration schemes 
in dissipative particle dynamics give rise to pronounced artifacts in physical quantities such as 
the compressibility and the diffusion coefficient. We assess the quality of these integration schemes, 
including variants based on a recently suggested self-consistent approach, and examine their relative 
performance. Implications of integrator-induced effects are discussed. 
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One of the current challenges in theoretical physics 
is to understand the basic principles that govern collec- 
tive properties of soft-matter systems. From a modeling 
point of view, these systems are problematic due to the 
fact that numerous phenomena take place at mesoscopic 
time and length scales, while the most accurate "brute- 
force" molecular dynamics simulations are limited to mi- 
croscopic time and length scales. To overcome this prob- 
lem, a number of "coarse-grained" approaches Jl|^ have 
been suggested and developed to simplify the underlying 
microscopic model while retaining the essential physics. 

Introduced in 1992 [|J and cast into its present form 
in 1995 Dissipative Particle Dynamics (DPD) has 
become one of the most promising methods for soft- 
matter simulations 00]- From a technical point of view, 
DPD differs from Molecular Dynamics (MD) in two re- 
spects. First, the conservative pairwise forces between 
DPD particles (which represent clusters of microscopic 
particles) are soft -repulsive, which makes it possible to 
extend the simulations to longer time scales. Second, a 
special "DPD thermostat" for the canonical ensemble is 
implemented in terms of dissipative as well as random 
pairwise forces such that the momentum is locally con- 
served, which results in the emergence of hydrodynamic 
flow effects on the macroscopic scale. 

However, the pairwise coupling of particles by the dis- 
sipative and random forces makes the integration of the 
equations of motion a non-trivial task. It has been ob- 
served that essentially all traditional integration schemes 
lead to distinct deviations from the true equilibrium be- 
havior, including an unphysical systematic drift of the 
temperature from the value predicted by the fluctuation- 
dissipation theorem, and artificial structures in the ra- 
dial distribution function fl|6|-|l0|. Consequently, vari- 
ous integration schemes have been suggested to overcome 
these problems Q^H- Some approaches are based on 
the use of phcnomcnological "tuning parameters" which 
mimic higher-order corrections in the integration pro- 
cedure BM . A more elaborate technique suggested by 



Pagonabarraga et al. || determines the velocities and 
velocity-dependent dissipative forces in a self-consistent 
fashion. Although both approaches have been shown to 
reduce numerical artifacts in some cases j|,[i].|]f| , there 
is still no good understanding as to which integration 
scheme is most suitable for future extensive soft-matter 
DPD simulations, and a thorough comparison including 
recently suggested schemes ||fj is pending. Moreover, 
the effect of integrators on dynamic quantities such as 
transport coefficients has received only little attention so 
far plj . In this work, we address these issues. 

We consider various integrators based on the velocity- 
Verlet scheme jl2| and assess their quality by studying 
a number of physical observables such as temperature, 
radial distribution function, compressibility, and tracer 
diffusion. We demonstrate that there is no reason to 
use integrators which contain tuning parameters, since 
better schemes are readily available. However, we also 
find that even a self-consistent approach gives rise to 
subtle temperature drifts, which can be corrected by a 
method presented in this work. 

For a system of N particles with mass to, coordinates 
{r^}, and velocities {v^}, the pairwise conservative, dis- 
sipative, and random forces exerted on particle "£" by 
particle "j" are given by, respectively, 
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The £ij are symmetric random variables with 
zero mean and unit variance, uncorrelated for different 
pairs of particles and different times. The forces are soft- 
repulsive due to the weight function w(rjj) for which we 



adopt the commonly made choice w(r^) = 1 
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j ■) = for r y > r c , with a cut-off 



distance r c |j4|. The strength of the conservative, dissipa- 
tive, and random forces is determined by the parameters 
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a, 7, and a, respectively. The equations of motion are 
then given by the set of stochastic differential equations 
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where Ff = YljjH ^fj ^ s * ne total conservative force 
acting on particle "z" (with Ff and Ff defined corre- 
spondingly). This continuous-time version of DPD sat- 
isfies detailed balance and describes the canonical ensem- 
ble if a and 7 obey the fluctuation-dissipation relation 
a 2 / 7 = 2fc B T*§. 

In order to integrate the equations of motion, we use 
the widely adopted velo city- Ver let scheme |l2j as a start- 
ing point and consider the most commonly used integra- 
tors based on this approach. These are summarized in 
Table Q, where the acronym "MD-VV" corresponds to 
the standard velocity- Verlet algorithm used in classical 
MD simulations. Unlike in MD, however, the forces in 
DPD depend on the velocities. For that reason Groot 
and Warren || proposed a modified velocity- Verlet in- 
tegrator ("GW(A)" in Table |). In this approach, the 
forces are still updated only once per integration step, 
but the dissipative forces are evaluated based on inter- 
mediate "predicted" velocities Vj . The calculation of 
involves the use of a phenomenological tuning parameter 
A which mimicks higher-order corrections in the integra- 
tion procedure. The problem is that the optimal value of 
A, which minimizes temperature drift and other artifacts, 
depends on model parameters and has to be determined 
empirically. Recently, Gibson et al. Q proposed a slightly 
modified version of the GW integrator. This "GCC(A)" 

TABLE I. Update schemes for a single integration step for 
various DPD integrators (acronyms see text). 

GW(A): steps (0)-(4),(s) 

MD-VV = GW(A = l/2) : steps (l)-(4), (s) a 

GCC(A): steps (0)-(5), (s) 

DPD-VV = GCC(A = l/2) : steps (l)-(5), (s) a 
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a with substitution of Vj for Vj in step (3). 

b Sampling step (calculation of temperature ksT , g(r), 



integrator updates the dissipative forces (step (5) in Ta- 
ble Q) for a second time at the end of each integration 
step. Choosing A =1/2 in the GCC integrator is equiva- 
lent to the MD-VV scheme supplemented by the second 
update of the dissipative forces. This Ver let-type in- 
tegrator, here termed "DPD-VV" , is appealing because 
it does not involve a tuning parameter, yet takes the 
velocity-dependence of the dissipative forces at least ap- 
proximately into account. 

Unfortunately, all of the above integrators display pro- 
nounced unphysical artifacts in g(r) and thus do not 
produce the correct equilibrium properties (see Fig. |l| 
and discussion below). This highlights the need for an 
approach in which the velocities and dissipative forces 
are determined in a self-consistent fashion. To this end, 
we present in Table || the update schemes for two self- 
consistent variants of DPD-VV. The basic variant, which 
is similar in spirit to the self-consistent leap-frog scheme 
introduced by Pagonabarraga et al. J8|, determines the 
velocities and dissipative forces self-consistently through 
functional iteration, and the convergence of the iteration 
process is monitored by the instantaneous temperature 
ksT. In the second approach, we furthermore couple 
the system to an auxiliary thermostat, thus obtaining an 
"extended-system" method in the spirit of Nose-Hoover 
fl3| (see below for details). 

Since the problems due to dissipative and stochastic 
forces in DPD are particularly noticeable in the absence 
of conservative forces, we therefore focus on a 3D 
ideal gas @,||. In our simulations we use a box of size 

TABLE II. Update scheme for self-consistent DPD-VV 
without and with (steps (i)-(iii)) auxiliary thermostat. The 
self-consistency loop is over steps (4b) and (5) as indicated. 
The desired temperature is ksT* . Initialization: r\ — 0, 
7 = <T 2 /(2fc.gT'*), and ksT calculated from the initial velocity 
distribution. 
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FIG. 1. Radial distribution functions g(r) as obtained in 
DPD simulations of the 3D ideal gas for At = 0.01 (left) and 
At — 0.1 (right) for velocity- Verlet-based integrators. Time 
is given in units of r c ^/m/kBT. 

10 x 10 x 10 with periodic boundary conditions, a ran- 
dom force strength a — 3 , and a particle density p = 4 
(i.e., N = 4000 particles) Q. For the equilibrated 
system, the temperature (k B T) and the radial distri- 
bution function g(r) were sampled. For the ideal gas, 
g(r) = 1 in the continuum limit, and therefore any devi- 
ation from 1 has to be interpreted as an artifact due to 
the employed integration scheme. Artifacts in g(r) are 
also reflected in the relative isothermal compressibility 
k t = n T /ntf ea \ where n^ eal = {pk B T*)- x denotes the 
compressibility of the ideal gas in the continuum limit. 
For an arbitrary fluid, Ht is related to g{r) by kt — 
1 + 4irp f^°dr r 2 [g(r) — 1] , and thus any deviation from 
kt = 1 for the ideal gas indicates an integrator-induced 
artifact. Finally, to gauge underlying problems in the ac- 
tual dynamics of the system, we consider the tracer diffu- 
sion coefficient D T = lim^oo Y,fLi([ r i( t ) ~ r i(°)] 2 )> 
which characterizes temporal correlations between the 
displacements (velocities) of the tagged particle. 

Results for g{r) are shown in Fig. [|. We find that the 
deviations from g{r) = 1 are very pronounced for MD- 
VV, indicating that even at small time steps it gives rise 
to unphysical correlations. The performance of DPD-VV 
is clearly better, while the self-consistent scheme leads 
to even smaller deviations. Studies of the integrators 
GW(A = 0.65) and GCC(A) (for a few values of A) re- 




FIG. 2. Left: Dimensionless isothermal compressibility kt 
vs. At for the integrators shown in the legend. Right: Tracer 
diffusion coefficient Dt vs. At for the same integrators. 
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FIG. 3. Double-logarithmic plot of the modulus of the de- 
viation of <fcsT> from the desired temperature ksT* = 1 vs. 
At. For self-consistent DPD-VV with auxiliary thermostat, 
la error bars are shown for some of the data points. 

vealed that their results were approximately similar to 
those of MD-VV and DPD-VV, respectively. For all 
integrators, the artificial structure in g(r) becomes more 
pronounced with increasing time increment At, and it is 
intriguing that the bias introduced by the self-consistent 
integrator for At = 0.10 is comparable to that introduced 
by MD-VV for At = 0.01. 

The relative isothermal compressibilities kt evaluated 
from g(r) are shown in Fig. |2[ The qualitative behavior 
of kt reflects our findings for g(r) Jig ]. However, the 
magnitude of deviations from kt = 1 is astounding, and 
raises serious concern for studies of response functions 
such as the compressibility for interacting fluids close to 
phase boundaries. Similarly, the results for tracer dif- 
fusion (also in Fig. ||) indicate that DPD-VV and the 
self-consistent approach work well up to reasonably large 
time steps, while the other integrators were found to per- 
form less well. Thus, the decay of velocity correlations in 
tracer diffusion is sensitive to the choice of the integrator. 
These results demonstrate that special care is needed in 
studies of DPD model systems, and suggest that integra- 
tors commonly used in MD should not be employed in 
DPD as such. 

Next we discuss the deviations of the observed actual 
temperature (k B T) from the desired temperature k B T* 
(see Fig. ||). For MD-VV this "temperature drift" is 
always positive and increases monotonically with At. For 
DPD-VV, (ksT) first decreases with increasing At, then 
exhibits a minimum at Atsa0.25, and eventually becomes 
larger than k B T* . The self-consistent approach exhibits 
a negative, monotonically increasing temperature drift 
up to At w 0.13, where this scheme becomes unstable 
at the employed particle density. Most importantly, we 
find, surprisingly, that the modulus of the temperature 
deviation is even larger than the one for DPD-VV. In 
a recent work, Pagonabarraga et al. || studied the 2D 
ideal gas using a self-consistent version of the leap-frog 
algorithm, and found good temperature control for At = 
0.06 at p — 0.5. This discrepancy can be explained by 
our observation for the 3D ideal gas that the temperature 
drift is in general more pronounced at higher densities. 
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In cases where temperature preservation is crucial 
in calculating equilibrium quantities, we finally demon- 
strate how this can be achieved. The idea is to supple- 
ment the self-consistent scheme by an auxiliary thermo- 
stat, which preserves the pairwise conservation of mo- 
mentum by employing a fluctuating dissipation strength 



7(*) = 



2k B T 



-(1+ 77(f) At), 



(3) 



where 77 is a thermostat variable. The rate of change of 77 
is proportional to the instantaneous temperature devia- 
tion, 77 = C '(ksT—kBT*) , where C is a coupling constant 
(step (i) in Table ||) (T§. This first-order differential 
equation must be integrated (step (ii)) simultaneously 
with the equations of motion. In this respect our ther- 
mostat resembles the Nose-Hoover thermostat familiar 
from MD simulations jl3| . 

For this extended-system method, we find (Fig. ||) that 
the temperature deviations diminish by over two orders 
of magnitude, with a modulus typically of the order of 
1CU 5 . . . 1CU 4 . We also found virtually the same results 
for g(r) and /?t as for the self-consistent scheme without 
the thermostat. This suggests that the auxiliary ther- 
mostat is useful in studies of equilibrium quantities such 
as the speficic heat. However, we feel that the auxiliary 
thermostat is not an ideal approach to describe quanti- 
tative aspects of tracer diffusion. A more detailed study 
is currently in progress Jig ]. 

In this work, we have shown that integration schemes 
may in DPD lead to pronounced artifacts in response 
functions and transport coefficients. This constitutes a 
serious problem for studies of soft systems, and high- 
lights the timely need to resolve this issue. We have 
demonstrated that these artifacts can be sufficiently sup- 
pressed by using velocity- Verlet-based schemes in which 
the velocity dependence of the dissipative forces is taken 
into account. The velo city- Ver let scheme without iter- 
ations but with an additional update of the dissipative 
forces (DPD-VV) performs — at essentially unchanged 
computational costs — already considerably better than 
the Groot-Warren integrator. The best overall perfor- 
mance is found for a recently proposed approach in 
which particle velocities and velocity-dependent dissipa- 
tive forces are determined self-consistently. This scheme 
provides an accurate description for the quantities stud- 
ied here, except for the temperature whose persisting 
drift is found to be unexpectedly significant. As shown 
in this work, however, this drift can be suppressed by at 
least two orders of magnitude by a scheme which supple- 
ments self-consistency with an auxiliary thermostat. The 
computational cost of the self-consistent scheme with- 
out thermostat remains modest, and increases only by a 
factor of 1.5 to 3 with respect to the standard velocity- 
Verlet algorithm. 

Although DPD has been very successful especially 
in simulations of polymeric systems and in reproduc- 
ing equilibrium properties, there have been doubts as to 



whether DPD is able to describe the dynamics and trans- 
port properties of complex fluids [|],|8| . Since the origin of 
the observed discrepancies is not clear and a general the- 
ory is still lacking, it will be interesting to see what the 
impact of the ideas presented here is on transport prop- 
erties. Work is in progress to address these questions 
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